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I. INTRODUCTION 



In bulk type-II superconductors, mixed states in which the vortices carry more than 
one flux quantum are unstable with respect to the usual vortices that carry a single flux 
quantum. Near the upper critical fleld this result was established by Abrikosov,- while near 
the lower critical fleld it is implied by Matricon's^ calculations for isolated vortices. As far as 
we are aware there have been no explicit calculations within Ginzburg-Landau (GL) theory 
for lattices of multiply quantized vortices at magnetic fleld values between these two limits; 
but since there is no reason to believe that such lattices should be stable at intermediate 
flelds there has been no reason to carry out the calculations. 

Our motivation for doing so comes from considering thin fllms of type-I superconductors. 
In a sufficiently weak perpendicular fleld the magnetic flux penetrates the sample in vortices 
rather than intermediate state structures, and in some circumstances multiply quantized 
vortices are present in Ginzburg-Landau theory calculations^'^ and experiments.^ We have 
been interested^ in working out the phase diagram for thin fllms of type-I GL supercon- 
ductors in arbitrary magnetic fleld, following Brandt's approach^*^ to accurate and efficient 
solutions of the GL equations. This approach involves iterative solutions of equations, and 
good initial values are essential because convergence is not guaranteed. It turns out that a 
solution of the GL equations in bulk can be used as starting point for the film problem — this 
is the case even for type-I superconductors, where the bulk vortex lattice solutions are unsta- 
ble with respect to the normal state. Hence, calculations for lattices of multiply quantized 
vortices in bulk superconductors are almost a prerequisite for the more physically motivated 
calculations in fllms. 

For brevity, in this paper "singles" and "doubles" are synonymous with lattices of singly 
and doubly quantized vortices (carrying one and two flux quanta, respectively). 

The formal developments we present are an extension of Brandt's work in Ref. [sj. In Sec. |TT] 
we describe how Brandt's Ansatz and iteration scheme for singles are modifled for doubles. 
In Sec. mil we offer illustrative results for order parameter and magnetic induction proflles, 
as well as the Gibbs free energy, for doubles and singles in a type-II superconductor. In 
the Appendix we discuss the solutions of the linearized GL equations, which serve as initial 
values for the iterative calculations. 
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II. FORMALISM AND IMPLEMENTATION 



Form of solutions for doubles 



Consider a vortex in which the phase of the order parameter changes by 27rp on circhng 
the vortex core. If that core is at the origin, then the order parameter behaves as ~ j-p^^vS 
as r — i- (see for example Tinkham^ Sec. 5.1 ) and the modulus squared order parameter 
as a; = ~ r^^. In this paper we focus on the p = 2 case. Brandt suggests^ that for a 
lattice with one vortex per primitive cell we adopt the Ansatz 
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[1 — COS (K ■ r)]' 



(1) 



in which r is a two-dimensional vector and K runs over reciprocal lattice vectors excluding 
the origin (as will be the case for all sums over K henceforth). This form satisfies the 
requirements of periodicity and fourth-power behavior near vortex cores. It turns out to be 
useful to express this with only first powers of cosines, 
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In a bulk superconductor B (r) = B (r) z. For small r the induction satisfies B (r) ^ 
B (0) — ^UJ (r), so B (0) — B (r) ~ r^. The small-r behavior suggests the following form for 
the deviation from mean induction, b{r) = B (r) — B, 



6(r) = 5^6K 
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2 cos (K ■ r) - - cos (2K • r) 



(3) 



The supervelocity can be decomposed as Q (r) = (r) + q(r), where Qa is the su- 
pervelocity in the Abrikosov limit, satisfying V x (r) = [i? — 2^q6 (r)] z, (with $o 
the flux quantum) and the deviation from the Abrikosov form has the property that 
z ■ (V X q (r)) = 6 (r). This last relation implies 
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zxK 



2 sin (K ■ r) - ^ sin (2K ■ r) 



(4) 



The mean induction B fixes the area S of the lattice unit cell, through 5* = 2^q/B. 
The unit cell has primitive lattice vectors Rio = XiX, Rqi = + and in those 
terms S = (Rio x Roi) ■ z = xi?/2- The general reciprocal lattice vector is K^n = 
2-71 [m?/2X + {mx2 + nxi) y] / S. 



B. Iterative expressions 



The GL free energy per unit volume F, referenced from the Meissner state, may be 
expressed in gauge invariant form and in the standard reduced units as^ 



2 



where angle brackets denote integration over a two-dimensional unit cell, (■ ■ ■ ) = 
^ fg ■ ■ - dx dy. Extremalization of F yields the first GL equation 

-V^cu = 2k^ [co -cu^ -uj\Qf - g] (6) 

where 

g = \ Vu\y4K^u (7) 

and the second GL equation 

V X Bz = -uQ (8) 

The first GL equation leads to an iterative equation for ok- Following Brandt, a stabi- 
lizing term 2k,'^u is added to both sides of Eq. ([6]) and then the equation is multiplied by 
cos K ■ r and integrated over the unit cell. The expansion ([2]) is inserted for u on the left side 
of the equation, and the orthogonality relation (cos (K ■ r) cos (K' ■ r)) = |5k,k' enables the 
integral on the left side to be done analytically. Rearranging leads to an identity which we 
treat step in an iterative solution for ok, 



Ok := 2 + u'^ + u |Q|^ + g) cosK ■ r) 

|K| +2«:2 



+ 4«K/2 

If K/2 is not a reciprocal lattice vector then aK/2 = 0; we will refer to such reciprocal lattice 
vectors as "fundamentals." Eq. (|9]) should be compared with the corresponding relation for 
singles, Eq. (11) in Ref. sj: the only differences are the existence of the second term and a 
factor of two in the first term. 

The next step in the iterative scheme is to rescale all of the so as to minimize F. This 
goes through without modification 

ttK ■■= a^iuj - u) \Qf - g) / {uj'^) (10) 



The iterative equation for 6k is derived in the same manner as Eq. ([9j). Taking the 
curl of Eq. ([8]), adding a stabihzing term —{u)B (r) to both sides, multiplying by cosK ■ r, 
integrating over the unit cell, inserting the expansion ([3]) on the left hand side, applying 
orthogonality of cosines, and rearranging leads to 

([(w - (a;))fi(r) + (Vw X Q) • z] cosK-r) 

+ 4^K/2 



|K|' + (a;) 



Again, if K is a fundamental then the coefficient 6k/2 = 0. 

The GL equations are solved, in principle, by cycling through Eqs. (I9l), ( fTOl) . and (fTTll 
until the coefficients converge to the desired level of precision. In practice we find that the 
equations as written do not usually converge to a physical solution; however, by "mixing" 
the Ok that comes out of ([9]) with the value from the prior iteration (and likewise for the 6k 
produced by f|TT]) ) the convergence of the algorithm is much improved, though at the cost of 
more iterations. We have not attempted to determine optimal mixing parameters. Taking 
90% of the prior iteration plus 10% of the current iteration is sufficient for every calculation 
we have carried out so far. 

Even with mixing, it is crucial to have a good initial guess for the ok and 6k- Brandt 
has demonstrated^ that the solution of the linearized GL equations for the ok, together 
with 6k = for all K, serves well for the initial values for singles, and we find the same to 
be true for doubles. Constructing solutions to the linearized GL equations for doubles in 
terms of the ok is not trivial, and we detail our method in Appendix |Al The solution to 
the linearized GL equations is used to begin the iteration cycle in Eq. (ITU]) . 



C. Implementation issues 

In order to have a finite computational problem the expansions for a;, 6, and q must 
be truncated; and the iterations involve integrals over the unit cell which must be numer- 
ically evaluated. These two issues are related. For singles, it is sufficient to carry out the 
quadrature by summation of values on a grid aligned with the primitive lattice vectors, and 
to include in the expansions only |K| < -ft'max with i^max chosen so that the number of 
reciprocal lattice vectors is the same as the number of points in the integration grid. For 
doubles the situation is more complicated. 
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Observe that the expansion 
form as for singles, 
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for u can be rearranged so that it has nearly the same 
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(12) 



When K is a fundamental, as defined following Eq. (jO]), aK/2 = 0. Eqs. ([T]), ([2]) and ffT2ll 
are identical for infinite sums, but they are different when truncated. As discussed in Ap- 
pendix |X1 the Ok that solve the linearized GL equations for doubles do not fall off in a 
Gaussian manner like they do for singles; however, the 2aK — |ctK/2 are approximately 
Gaussian. This motivates the following truncation scheme for constructing u when eval- 
uating the integrals over the unit cell: use Eq. [12], including in the sum reciprocal lattice 
vectors with |K| < -ft"max except for fundamentals with i^max/2 < |K| < Kmax- Expressions 
analogous to Eq. [T2] exist for b and q, namely 
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(14) 



and we apply the same truncation scheme. 

We conclude this section with a warning. It is tempting to define Ck = 2aK — 
d-K = 26k — \bvi/2) and carry out iterative calculations for those quantities, by moving the 
\0''K/2 from the right side of Eq. ([9]) to the left (and likewise for Eq. f|TT]) ): the resulting 
equations have exactly the form of Eqs. 11 and 13 from Ref. js] for the iterations of the 
coefficients for singles. Doing this invariably leads to a "singles-like" solution {u ~ and 
B (0) — B (r) ~ for small r) which is inconsistent with the assumed forms of B, Qa and 
S; in addition this unphysical solution is a free energy saddle point rather than a minimum. 



III. RESULTS AND CONCLUSIONS 



To illustrate the effectiveness of the calculational scheme described in the previous section, 
we will present results for k, = 1 for triangular arrays of both doubly and singly quantized 
vortices. In comparing doubles with singles it is important to do the calculations on an 
equal footing; with this in mind, in order to have the same spacing between real-space 
grid points in both calculations we use about twice as many grid points for doubles than 
for singles. For B greater than hqHc2/5 we typically use 32 points along each primitive 
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lattice vector for singles and 46 for doubles. Singles and doubles then have the same -ft'max 
but the doubles calculation includes twice as many reciprocal lattice vectors as the singles 
calculation. A typical singles calculation converges in approximately 50 iterations while 
the doubles require about four times as many. At lower inductions the area of the vortex 
core becomes considerably less than the area of the unit cell, and in order to represent the 
solutions well the real-space sampling needs to be refined by increasing the number of grid 
points per unit cell edge to 64 or 96 for singles and with proportionally increased numbers 
for doubles. The growth in the number of grid points and reciprocal lattice vectors puts a 
practical lower limit on the mean induction of /io-f^c2/10. 

In Figure [T] we show the order parameter and induction along a line connecting two 
adjacent vortices a.t B = /ioi/c2/10 for both doubles and singles. As one would expect, the 
vortex cores for doubles are wider than for singles. In Fig.[2]we present the Gibbs free energy 
density G = F — ■ B in the full range of applied fields Hd < Ha < i^c2 for triangular 
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singles and doubles. The applied field Ha is calculated in the same manner as in Refs. l8 
based on the virial theorem of Doria, Gubernatis, and Rainer.— These results confirm that 
doubles are thermodynamically unstable in bulk type-II superconductors. Similar results 
are obtained for other values of /t > 1/ \/2. The calculations also yield convergent results for 
type-I superconductors but in such cases all vortex states are unstable with respect to the 
Meissner state. 

In conclusion, we have produced precise numerical solutions to the GL equations consist- 
ing of infinite lattices of "doubly quantized" vortices in bulk superconductors. The calcula- 
tions can be carried out efficiently for mean inductions down to 10% of the upper critical 
value. Although such solutions of the GL equations never globally minimize the GL free 
energy for bulk superconductors, we expect they will be useful as starting points for solving 
the GL equations in film geometry. 



Appendix A: Solving the Linearized GL equations in Terms of the ok 

In his pioneering work on vortex lattices in superconductors, Abrikosovi showed that for 
an applied field just below ifc2, the first GL equation (when expressed as an equation for 
the order parameter) has the form of Schrodinger's equation for a charged particle confined 
to a plane and subject to a magnetic field. With an assumed periodicity of the vortex lines 
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and one flux quantum per vortex, an analytic solution ipA exists and can be expresseciiii^iii^ 
in terms of a Jacobi theta function, 

{x, y) = e-^y'l^-vH^ {x + ly) , ^-1±J^ (Al) 

where {z, r) = 2 Yl'^=o (^)" e*'^^^"''"^) sin(2n + l)z and the lattice parameters Xi, X2 and 
y2 were deflned just below Eq. With this form for ipA its modulus squared, Ua = I'^pA^, 
is expressed as a double sum. This leads^ to the Fourier like expansion of real terms, 
CO A = YIk [-'- ~ "^^^ (-^ ■ ^)]' 'with Ok ~ ^ i^_-^m+mn+n ^-K^^s/8n ^ Notc that the sum over 
K is still a double sum over m and n. 

Lasher— pointed out that for vortices of multiplicity p, tjj^^^ (r) = ^ipA {^/^/p)Y ^ 
corresponding solution of the linearized GL equations. In principle one could use this form 
to determine for doubles in the expansion ([1]), starting from (lAip . but we did not attempt 
to carry that through. 

We have taken an alternative approach based on numerical solution of a linear system for 
the Ok derived from the linearized GL equations. In the linear regime b (r) = -B(O) — B — 
loa (r) /2k (see, for example, De Gennes^ Sec. 6.7). Taking the curl and combining with the 
second GL equation yields 

■^VujA X z = uaQa (A2) 

Combining (1A2I) , ([2]) , and the Fourier expansion for the supervelocity in the Abrikosov limit 
(see Eq. (24) in Ref. 3), then applying uniqueness of Fourier series leads to the linear system 

Y,A,,ai^ = -{uA) (A3) 

with 
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We follow Brandt's convention that {uja} = 1, so (3/2) '^k ~ ^■ 

The infinite system of equations (1A3P could be rendered finite by setting = for 
|K| > Kmax; however, this is not a good closure assumption because of slow convergence 
with increasing -ft'max- Eq. (lA4p shows the strong connection between and 0^/2 Previously 
mentioned in Sec. Ill C[ In particular, for < |K| < i^inax the corresponding 

are connected to coefficients associated with vectors beyond the cutoff. In addition, as 
|Kj| — 7- 00, Cji —7- —3/2, which leads to ^ large |K|. (For large |Kj|, 

YIk ^ji'^Ki ~ ~(3/2)X]k'^k ~ ~{<^a), hence the latter two terms on the right side of 
Eq. IA4I must sum to zero.) We therefore set a^j^ = a^/4, afj^ = a^/16, and so on for 
-^max/2 < |K| < Kmax and this leads to a modified linear system with coefficients A^j. For 
i such that < |Kj| < /Tmax, 

00 

A;, = A,, + 5^4-'C,.2., (A6) 
1=1 

with K2ij = 2'Kj. We truncate the sum at / = 4 after finding no change in the results for 
when further terms are included. 

In Figs. |3] and H] we show the results of numerically solving the linear system for k = 1 
and -ft'max = 36. In the former only for fundamental reciprocal lattice vectors are shown; 
while in the latter for reciprocal lattice vectors that are powers of two times several 
different fundamentals are displayed, showing that ~ holds even for |K| not very 

large. 

Finally, let us note that as a check on this approach to solving the linearized GL equations 
we have carried an analogous analysis for singles. The numerical results from solving the 
corresponding linear equation for the match the exact, analytic results. 
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Figure 1: Cross sections along y = of induction (curves decreasing from x = 0) and squared 
order parameter (curves increasing from x = 0) for doubly quantized (dashed) and singly quantized 
(solid) triangular vortex lattices at the same parameters k = 1 and B = /xo-ffc2/10. The inter-vortex 
spacing for singles is 8. 5 A and is greater by a factor of \/2 for doubles. 
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Figure 2: Gibbs free energy density, referenced from the normal state, as a function of applied field 
at GL parameter k = 1. 
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Figure 3: Solutions to the linear system for the associated with fundamental reciprocal lattice 
vectors. 
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Figure 4: Solutions to the linear system for the ag^ associated with eight fundamental reciprocal 
lattice vectors and those vectors multiplied by powers of two. The lines are guides to the eye 
connecting a^, o^K' ^4k'- • • 
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